---
title: "05 Propensity Score Weighting Data-Antenatal"
author: "Ye Shen, Radhika Tampi, Raj Vatsa"
date: "11/27/2021"
output: pdf_document
---
Balancing weights were created using household wealth index, a time-varying confounding variable that substantially differed between intervention and control in 2014 and 2016. 
```{r, results = 'hide', message = FALSE}
knitr::opts_chunk$set("warning"=FALSE)
library(tidyverse)
library(survey)
library(WeightIt)
library(survey)
library(questionr)#use survey design in ggplot
library(patchwork)
library(cobalt)
library(gt)
library(optweight)
library(mlogit)
```
## Read Cleaned Data and Prepare for Balance Weighting
Create a categorical variable of group+year combo for each dataset to apply multinomial weighting using method proposed by Stuart et al 2014 for integrating DiD and weighting
```{r}
antenatal_did<-readRDS("data/cleaned_antenatal_did.rds")%>%
                          mutate(group_year=as.factor(case_when(
                            group==1&year==0 ~ 1, #Treatment in pre-period: Group 1 in Stuart et al 2014
                            group==1&year==1 ~ 2, #Treatment in post-period: Group 2 in Stuart et al 2014
                            group==0&year==0 ~ 3, #Control in pre-period: Group 3 in Stuart et al 2014
                            group==0&year==1 ~ 4 #Control in post-period: Group 4 in Stuart et al 2014
                          )))

 
 svy_antenatal_did <- svydesign(
    id = ~cluster, strata = ~group, weights = ~wmweight,
    data = antenatal_did
  )
```
## Before balance weighting
```{r}
svy_antenatal_14<-subset(svy_antenatal_did, year==0)
svy_antenatal_16<-subset(svy_antenatal_did, year==1)
```

### Distributions of wealth index score between treatment and control in 2014 and 2016
```{r, message=FALSE}
antenatal_wealth_2014plot<-ggsurvey(design = svy_antenatal_14, mapping = aes(x=as.factor(group) , y=menage_score_bien_etre))+
  geom_violin()+
  scale_y_continuous(name="Wealth Index Score", limits=c(-1,5),breaks=seq(-1,5,1))+
  scale_x_discrete(name="", labels=c("0"="Control", "1"="Treatment"))+
  theme_bw()+
  labs(title="2014")
antenatal_wealth_2016plot<-ggsurvey(design = svy_antenatal_16, mapping = aes(x=as.factor(group) , y=menage_score_bien_etre))+
  geom_violin()+
   scale_y_continuous(name="Wealth Index Score", limits=c(-1,5),breaks=seq(-1,5,1))+
  scale_x_discrete(name="", labels=c("0"="Control", "1"="Treatment"))+
  theme_bw()+
  labs(title="2016")
antenatal_wealth_2014plot+antenatal_wealth_2016plot +plot_annotation(tag_levels=c("A"),title="Fig X1. Distributions of wealth index score between treatment and control \n in pre and post- periods among pregnant women")
```
### Household Wealth Index Imbalance QQPlot 
```{r, message=FALSE}
gg_qqplotimbalance<-function(a, b, quantiles = seq(0, 1, 0.01)){

  
  ggplot(mapping = aes(x = quantile(a, quantiles), 
                       y = quantile(b, quantiles))) + 
    geom_point() +
    geom_abline(aes(slope = 1, intercept = 0), linetype = 2) 
    
}

# svy_antenatal_16_trt<-subset(svy_antenatal_16, group==1)
# svy_antenatal_16_con<-subset(svy_antenatal_16, group==0)

antenatal_wealth_2016qqplot<-gg_qqplotimbalance(subset(svy_antenatal_16, group==1)$variables$menage_score_bien_etre,subset(svy_antenatal_16, group==0)$variables$menage_score_bien_etre)+
      labs(x = "Treatment Wealth Index Score",
         y = "Control Wealth Index Score",
         title="2016")+
          scale_y_continuous(limits = c(-1,5))+
          scale_x_continuous(limits = c(-1,5))+
          theme_classic()
          
                              
antenatal_wealth_2014qqplot<-gg_qqplotimbalance(subset(svy_antenatal_14, group==1)$variables$menage_score_bien_etre,subset(svy_antenatal_14, group==0)$variables$menage_score_bien_etre)+
      labs(x = "Treatment Wealth Index Score",
         y = "Control Wealth Index Score",
         title="2014")+
          scale_y_continuous(limits = c(-1,5))+
          scale_x_continuous(limits = c(-1,5))+
          theme_classic()

antenatal_wealth_2014qqplot+antenatal_wealth_2016qqplot +plot_annotation(tag_levels=c("A"),title="Fig X2. Imbalance qqplot on wealth index score between treatment and control \n in pre and post- periods among pregnant women")
```

### Combine the two imbalance plots into one figure
```{r}
(antenatal_wealth_2014plot+antenatal_wealth_2016plot) / (antenatal_wealth_2014qqplot+antenatal_wealth_2016qqplot)+plot_annotation(tag_levels=c("A"),title="Fig S1. Imbalance in household wealth index score between treatment and control \n in pre and post- periods among pregnant women", caption="Note: Panels A & B are violin plots displaying the wealth index score distributions for each group. \n Panels C & D are qq plots contrasting distributions in treatment against control.")&
    theme(plot.tag = element_text(size = 10), text = element_text( face="bold"))

ggsave("table_fig/paper/Supplemental/fig_s1_imbal_wealth_antenatal.png",width = 200,height=200, units="mm")
```

## Balance Weighting

### Optimization-Based Weights
Method with best performance re: smallest coefficients of variance close to 0 and large effective sample sizes. We used this method to create weights for the extension analysis, but showed how the second best method and the defaul method performed below.
```{r}
Weighted_antenatal_opt<-weightit(group_year~menage_score_bien_etre,data=antenatal_did, s.weights = antenatal_did$wmweight, method="optweight")
summary(Weighted_antenatal_opt)


cobalt::bal.tab(Weighted_antenatal_opt)


bal.plot(Weighted_antenatal_opt,"menage_score_bien_etre", which="both")+
  labs(title="Fig S5. Balance among 4 treatment-by-year groups before and after weighting on \n wealth index score among pregnant women", x="Wealth Index Score", fill="Treatment and Period")+
  scale_fill_discrete(labels=c("1"="Treatment Preperiod","2"="Treatment Postperiod","3"="Control Preperiod","4"="Control Postperiod"))+
  theme(plot.title = element_text(size=14, face="bold"))
ggsave("table_fig/paper/Supplemental/fig_s5_weighting_balance_antenatal.png")
```
The (standardized) difference in means between the two groups after adjusting using Optimization-Based Weights is zero.

### Entropy Balancing 
Close second best performance re: smallest coefficients of variance close to 0 and large effective sample sizes

Note: This method seems to be the same as "ebcw" Empirical Balancing Calibration Weights
```{r}
Weighted_antenatal_ebal<-weightit(group_year~menage_score_bien_etre,data=antenatal_did, s.weights = antenatal_did$wmweight, method="ebal")
summary(Weighted_antenatal_ebal)
cobalt::bal.tab(Weighted_antenatal_ebal)


bal.plot(Weighted_antenatal_ebal,"menage_score_bien_etre", which="both")+
  labs(title="FigX4. Balance among 4 treatment-by-year groups before and after weighting \n on wealth index score for pregnant women")+
  scale_fill_discrete(labels=c("1"="Treatment Preperiod","2"="Treatment Postperiod","3"="Control Preperiod","4"="Control Postperiod"))
```
### Multinomial regression PS (default)
```{r}
Weighted_antenatal_ps<-weightit(group_year~menage_score_bien_etre,data=antenatal_did, s.weights = antenatal_did$wmweight, method="ps")
summary(Weighted_antenatal_ps)
cobalt::bal.tab(Weighted_antenatal_ps)


bal.plot(Weighted_antenatal_ps,"menage_score_bien_etre", which="both")+
  labs(title="FigX4. Balance among 4 treatment-by-year groups before and after weighting \n on wealth index score for pregnant women", x="Wealth Index Score")+
  scale_fill_discrete(labels=c("1"="Treatment Preperiod","2"="Treatment Postperiod","3"="Control Preperiod","4"="Control Postperiod"))
```



## After balance weighting using optweight
### Assign survey design to did dataset using the new weights
```{r}
antenatal_did$balwt<-Weighted_antenatal_opt$weights

svy_antenatal_did_wted <- svydesign(
    id = ~cluster, strata = ~group, weights = ~balwt,
    data = antenatal_did
  )
```
### Subset to each year for weighted summary stats in table 2
```{r}
#subset did survey data to 2014 and 2016 for summary stats
  svy_antenatal_2014 <- subset(svy_antenatal_did_wted,year==0)
  svy_antenatal_2016 <- subset(svy_antenatal_did_wted,year==1)
```

### Weighted DiD regression adjusting for wealth index score; adjusted and unweighted; unajusted and unweighted
 * Outcome: Any antenatal consultation
```{r}
#Adjusted and weighted DID regression on Any antenatal consultation
  did_antenatal_wtd <- svyglm(antenatal.publichealthcenter.consultation~year*group+menage_score_bien_etre, design = svy_antenatal_did_wted, family = gaussian)

#Unadjusted and weighted DID regression on Any antenatal consultation
  did_antenatal_wtd_unadj <- svyglm(antenatal.publichealthcenter.consultation~year*group, design = svy_antenatal_did_wted, family = gaussian)
  
#Unadjusted and unweighted DID regression on Any antenatal consultation
  did_antenatal_unwtd_unadj <- svyglm(antenatal.publichealthcenter.consultation~year*group, design = svy_antenatal_did, family = gaussian)
```

## Table for pregnant women
 * Outcome: Any antenatal consultation
### Prepare Table
 

```{r}
table3_antenatal <- 
  as.data.frame(array(data = NA, dim = c(2, 22)))
# "Indicator","Study group", "Denominator", "Numerator", "Percent (SE)", "Difference 
# (p-value)","Denominator", "Numerator", "Percent (SE)", "Difference (p-value)","Percent 
# relative change", "Percent absolute change", "Difference in differences (p-value)"
colnames(table3_antenatal) <- 
  c("indicator", "group", "d_2014", "n_2014", "m_2014", "se_2014", "dif_2014", "p-val_2014", 
    "d_2016", "n_2016", "m_2016", "se_2016", "dif_2016", "p-val_2016", "trends_rel", 
    "trends_abs", "dif.in.dif(coef)_unadj_unwt", "p-val_did_unadj_unwt","dif.in.dif(coef)_unadj_wt", "p-val_did_unadj_wt","dif.in.dif(coef)_adj_wt", "p-val_did_adj_wt")
table3_antenatal$indicator <- "antenatal.publichealthcenter.consultation"

# Any visit (Access)
## Summary stats for 2014
table3_antenatal[rev(c(1,2)), 2] <- c(0,1)
table3_antenatal[rev(c(1,2)), 3] <- 
    round(svytable( ~group, design = svy_antenatal_2014)) # 'd_2014'
table3_antenatal[rev(c(1,2)), 4] <- 
    round(coef(svyby(~antenatal.publichealthcenter.consultation, by = ~group, design = svy_antenatal_2014, FUN = svytotal, keep.var = T, na.rm = T))) # 'n_2014'
table3_antenatal[rev(c(1,2)), 5] <- 
    round(coef(svyby(~antenatal.publichealthcenter.consultation, by = ~group, design = svy_antenatal_2014, FUN = svymean,keep.var = T, na.rm = T))*100,2) # 'm_2014'
table3_antenatal[rev(c(1,2)), 6] <- 
    round(SE(svyby(~antenatal.publichealthcenter.consultation, by = ~group, design = svy_antenatal_2014, FUN = svymean,keep.var = T, na.rm = T))*100,2) # 'se_2014'
table3_antenatal[rev(c(1,2)), 7] <- table3_antenatal[1, 5] - table3_antenatal[2, 5]# 'dif_2014'

### Estimate p-value of the 2014 difference with a Chi2 test
  if (table3_antenatal[1, 4] == 0 & 
      table3_antenatal[2, 4] == 0) {
    table3_antenatal[2, 8] <- NA
  } else {
    table3_antenatal[2, 8] <- 
      round(svychisq(~antenatal.publichealthcenter.consultation+group, design = svy_antenatal_2014)$p.value,2)
  }

## Summary stats for 2016
table3_antenatal[rev(c(1,2)), 9] <- 
    round(svytable( ~group, design = svy_antenatal_2016)) # 'd_2016'
table3_antenatal[rev(c(1,2)), 10] <- 
    round(coef(svyby(~antenatal.publichealthcenter.consultation, by = ~group, design = svy_antenatal_2016, FUN = svytotal,keep.var = T, na.rm = T))) # 'n_2016'
table3_antenatal[rev(c(1,2)), 11] <- 
    round(coef(svyby(~antenatal.publichealthcenter.consultation, by = ~group, design = svy_antenatal_2016, FUN = svymean, keep.var = T, na.rm = T))*100,2) # 'm_2016'
table3_antenatal[rev(c(1,2)), 12] <- 
    round(SE(svyby(~antenatal.publichealthcenter.consultation, by = ~group, design = svy_antenatal_2016, FUN = svymean,keep.var = T, na.rm = T))*100,2) # 'se_2016'
table3_antenatal[rev(c(1,2)), 13] <- table3_antenatal[1, 11] - table3_antenatal[2, 11]# 'dif_2016'

### Estimate p-value of the 2016 difference with a Chi2 test
  if (table3_antenatal[1, 10] == 0 & 
      table3_antenatal[2, 10] == 0) {
    table3_antenatal[2, 14] <- NA
  } else {
    table3_antenatal[2, 14] <- 
      round(svychisq(~antenatal.publichealthcenter.consultation+group, design = svy_antenatal_2016)$p.value,2)
  }

## Estimate values for trends 2014-2016
  table3_antenatal[rev(c(1, 2)), 15] <- 
    round((table3_antenatal[rev(c(1, 2)), 11] - 
       table3_antenatal[rev(c(1, 2)), 5]) / 
    table3_antenatal[rev(c(1, 2)), 5]*100,2) # 'trends_rel'
  table3_antenatal[rev(c(1, 2)), 16] <- 
    table3_antenatal[rev(c(1, 2)), 11] - 
    table3_antenatal[rev(c(1, 2)), 5] # 'trends_abs'

## DiD stats 
  #Unadj without balancing weight on wealth 
    #coef
  table3_antenatal[2, 17] <- round(summary(did_antenatal_unwtd_unadj)$coefficients[4, 1]*100,2)
    #P
  table3_antenatal[2, 18] <- round(summary(did_antenatal_unwtd_unadj)$coefficients[4, 4],2)
  #Unadj with balancing weight on wealth
    #coef
  table3_antenatal[2, 19] <- round(summary(did_antenatal_wtd_unadj)$coefficients[4, 1]*100,2)
    #P
  table3_antenatal[2, 20] <- round(summary(did_antenatal_wtd_unadj)$coefficients[4, 4],2)
  #Adj with balancing weight on wealth
    #coef
  table3_antenatal[2, 21] <- round(summary(did_antenatal_wtd)$coefficients[5, 1]*100,2)
    #P
  table3_antenatal[2, 22] <- round(summary(did_antenatal_wtd)$coefficients[5, 4],2)
```
### Final Supplemental Table 5 Antenatal
```{r}
# Replicate final Table 2
final_table3_antenatal<-data.frame(group=table3_antenatal$group, d_2014=table3_antenatal$d_2014, d_2016=table3_antenatal$d_2016,trends_rel= table3_antenatal$trends_rel,trends_abs= table3_antenatal$trends_abs)
final_table3_antenatal$m_2014 <- sprintf("%.2f (%.2f)", table3_antenatal$m_2014, 
                                       table3_antenatal$se_2014)
final_table3_antenatal$m_2016 <- sprintf("%.2f (%.2f)", table3_antenatal$m_2016, 
                                       table3_antenatal$se_2016)
final_table3_antenatal$dif_2014 <- ifelse(is.na(table3_antenatal$`p-val_2014`),"--" ,sprintf("%.2f (%.2f)", table3_antenatal$dif_2014,table3_antenatal$`p-val_2014`))
final_table3_antenatal$dif_2016 <- ifelse(is.na(table3_antenatal$`p-val_2016`),"--" ,sprintf("%.2f (%.2f)", table3_antenatal$dif_2016,table3_antenatal$`p-val_2016`))



final_table3_antenatal$`dif.in.dif(unadj_unwt)` <- ifelse(is.na(table3_antenatal$`p-val_did_unadj_unwt`),"--" , 
                                                        sprintf("%.2f (%.2f)", 
                                                                  table3_antenatal$`dif.in.dif(coef)_unadj_unwt`,
                                                                table3_antenatal$`p-val_did_unadj_unwt`))
final_table3_antenatal$`dif.in.dif(unadj_wt)` <- ifelse(is.na(table3_antenatal$`p-val_did_unadj_wt`),"--" , 
                                                        sprintf("%.2f (%.2f)", 
                                                                table3_antenatal$`dif.in.dif(coef)_unadj_wt`, 
                                                                table3_antenatal$`p-val_did_unadj_wt`))
final_table3_antenatal$`dif.in.dif(adj_wt)` <- ifelse(is.na(table3_antenatal$`p-val_did_adj_wt`),"--" ,
                                                sprintf("%.2f (%.2f)", 
                                                  table3_antenatal$`dif.in.dif(coef)_adj_wt`, 
                                                  table3_antenatal$`p-val_did_adj_wt`))


final_table3_antenatal <- as_tibble(final_table3_antenatal) %>% 
  mutate(
    group = case_when(group == 1 ~ "IG", 
                      group == 0 ~ "NG")
   ) %>%
  rename("orig_DiD" = "dif.in.dif(unadj_unwt)","unadj_wt_DiD" = "dif.in.dif(unadj_wt)", "adj_wt_DiD" = "dif.in.dif(adj_wt)")
  
gt(final_table3_antenatal) %>%
  tab_header(title = "Table S5 (Balance Weighting). Summary statistics of responses and trends on any antenatal visits at a PHF from the 2014 and 2016 IHOPE surveys") %>%
  tab_style(style = cell_text(align = "left"), location = cells_title("title")) %>%
  tab_spanner(
    label = "2014",
    columns = c(d_2014, m_2014, dif_2014)
  ) %>%
  tab_spanner(
    label = "2016",
    columns = c(d_2016, m_2016, dif_2016)
  ) %>%
  tab_spanner(
    label = "2014-2016 Trend",
    columns = c(trends_rel, trends_abs, orig_DiD,unadj_wt_DiD,adj_wt_DiD)
  ) %>%
  tab_style(style = cell_text(weight = "bold"), 
            location = cells_column_spanners(spanners = everything())) %>%
    cols_label(group = "Study group", d_2014 = "N", m_2014 = "Percent (SE)", 
               dif_2014 = " Difference (p-value)", d_2016 = "N", 
               m_2016 = "Percent (SE)", dif_2016 = " Difference (p-value)", 
               trends_rel = "Percent relative change", 
               trends_abs = "Percent absolute change", orig_DiD="Original Unadjusted DiD (p value) ",unadj_wt_DiD="Unadjusted DiD with balance weighting (p value)", adj_wt_DiD="Adjusted DiD with balance weighting (p value)") %>%
    tab_style(style = cell_text(weight = "bold"), 
              location = cells_column_labels(columns = everything())) %>%
    gt::gtsave("table_fig/paper/Supplemental/table_s5_antenatal_wt.png")
```

## Figure Antenatal
```{r}
num_var <- 1

fig_antenatal_wt <- as.data.frame(matrix(nrow = num_var * 4, ncol = 6))

fig_antenatal_wt$V1 <- round(c(
  table3_antenatal$m_2014, table3_antenatal$m_2016), 2)
fig_antenatal_wt$V2 <- round(c(
  table3_antenatal$n_2014, table3_antenatal$n_2016), 0)
fig_antenatal_wt$V3 <- round(c(
  table3_antenatal$d_2014, table3_antenatal$d_2016), 0)
fig_antenatal_wt$V4 <- rep(c("Intervention", "Non-Intervention"), num_var * 2)
fig_antenatal_wt$V5 <- rep(c("2014 survey", "2016 survey"), each = num_var * 2)
fig_antenatal_wt$V6 <- rep("Any Antenatal Visits to a PHF or CHW", num_var * 4) 
colnames(fig_antenatal_wt) <- c("mean", "num", "dem", "group", "year", "indicator")
fig_antenatal_wt$group <- factor(fig_antenatal_wt$group, levels = c("Non-Intervention", "Intervention"))
fig_antenatal_wt <- 
  ggplot(fig_antenatal_wt, 
         aes(x = mean, y = group, group = interaction(group, indicator), 
             label = sprintf("%s/%s", num, dem))) +
  geom_point(aes(colour = year), size = 3.5) +
  scale_colour_manual(breaks = fig_antenatal_wt$year, 
                      values = c("2014 survey" = "black", "2016 survey" = "red")) +
  geom_line() +
  labs(
    x = "Percentage", y = "",
    title = "Unadjusted with Balancing Weight on Wealth: \n Any Antenatal Visits to a PHF or CHW"
  ) +
  theme_bw() +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  scale_y_discrete(position = "right") +
  theme(
    legend.position = "bottom", axis.text = element_text(size = 10), 
    panel.grid.minor = element_blank(), legend.title = element_blank(),
    title = element_text(size = 12, face = "bold"), 
    axis.title = element_text(size = 13),
    legend.text = element_text(size = 12)
  ) + 
  annotate("text", label="{", x=50, y=1.5, size =20, fontface="plain")+
  annotate("text", label=paste0("Adjusted DiD (p-val): \n",final_table3_antenatal$adj_wt_DiD[2]), x=35, y=1.5)

saveRDS(fig_antenatal_wt,"table_fig/antenatal/fig_ante_wt.rds")
ggsave("table_fig/antenatal/fig3antenatal_wt.png")
```
